NON-EQUILIBRIUM THEORY OF THE ALLELE FREQUENCY 

SPECTRUM 

O ' STEVEN N. EVANS i, YELENA SHVETS 2, AND MONTGOMERY SLATKIN ■ 

O 

(N 



Abstract. A forward diffusion equation describing the evolution of the allele 
frequency spectrum is presented. The influx of mutations is accounted for 
by imposing a suitable boundary condition. For a Wright-Fisher diffusion 
with or without selection and varying population size, the boundary condition 

Vj . is liiUxiQ xf{x,t) = 9p{t), where f(-,t) is the frequency spectrum of derived 

alleles at independent loci at time t and p{t) is the relative population size 
at time t. When population size and selection intensity are independent of 
time, the forward equation is equivalent to the backwards diffusion usually 

^H ■ used to derive the frequency spectrum, but this approach allows computation 

of the time dependence of the spectrum both before an equilibrium is attained 
and when population size and selection intensity vary with time. From the 

,,^J diffusion equation, a set of ordinary differential equations for the moments of 

/(■,t) is derived and the expected spectrum of a finite sample is expressed in 
terms of those moments. The use of the forward equation is illustrated by 
considering neutral and selected alleles in a highly simplified model of human 

*vj \ history. For example, it is shown that approximately 30% of the expected total 

►^ , heterozygosity of neutral loci is attributable to mutations that arose since the 

^•^ • onset of population growth in roughly the last 150, 000 years. 
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ine allele- irequency spectrum is the distribution oi allele ire quencies at a large 
numb er of equivalent loci. The term "site-frequency spectrum" iJBraverman et al.l 

r^ , 119951) . is equivalent but emphasizes the application to individual nucleotides rather 

than alleles at different genetic loci. Here, we will use the frequency spectrum for 
both terms. 

^ ' Although models assuming reversible mutation predict an equilibrium distri- 

bution of allele frequencies (; Wright . 1931) . all recent studies of frequency spec- 

;J] ' tra assume irreversible mutation. Under that assumption, an equilibrium is not 

5^ \ reached at any locus, but the distribution across polymorphic loci reaches an equi- 

librium if both population size and selection intensities are constant. The theory 
predicting th e frequency spe ctru m under irrever si ble muta t ion \y as developed by 
iFisheJ l|l930() . IWrigyi)l938|) . and lKimural l|l964 - iKimural l|l969|) noted that this 
theory was a pplicable to nucleot i de po sitions and introduced the "infinitely many 
sites model." ISawver and Hartj l|l992l) incorporated the theory of Fisher, Wright 
and Kimura into a Poisson random field model for the purpose of estimating the 
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selection intensity from the observed frequency spectrum in a finite sample of chro- 
mos omes. Their method has been tested and refined bv iB^istamante et al.1 (|200ll) 
andlWilliamso n et all l)2004J) . 

Past population growth affects the frequency spectrum. iNei et al.l 1|1975(1 showed 
that ra pid g:ro\Y t h resu lted in more low frequency alleles than expected under neu- 
trality. iTaiimal l|l98 94 confirmed that conclusion and examined the eff ect of past 
popul ation growth on other aspects of the frequency spectrum. Griffiths and Tavara 
l)1998(l developed the coalescent theory for the frequency spectrum of neutr al alleles 
in a population that has experienced arbitrary changes in population size. Nielsen! 
112000') implemented the iGriffiths and Tavare (f998) simulation method and applied 
it to human SNP dat a for the purpo se of estimating the growth rate of human pop- 
ulations. Although iNielsenI 1)200(1 1 was not able to reject the hvpothesis of no 
growth, he noted that his analysis was of a small data set. IWakelev et al.l (|200ll) 
considered the same problem and develope d a method that allow ed for both ascer- 
tainment bias and population subdivision. Wakelev et al. (200f ) analyzed a larger 
SN P data set and found eyidenc e of recent growth. Wooding and Rogers (2002() 
and IPolanski and KimmeJ pOOiJ) also modeled the coalescent process underlying 
the spectrum of neutral alleles and developed analytic theory that allows for exact 
ca lculation of the spectrum for large sample sizes. 

iGriffitha l|2003j) reviewed and extended the theory of the frequency spectrum, 
making clear the role of time-reversal. He generalized that theory in two ways. He 
showed that the spectrum in a finite sample could be obtained from the solution 
to a backwards equation by assuming sampling with replacement, and he showed 
that the frequency spectrum in a population of variable size could be derived from 
the spectrum for a population of constant size when a transformation of the time 
scale reduces the backwards equation to one for a population of constant size. 
The transformation of time scales is always possible for neutral alleles, in which 
case the freguencv spectru m in a finite population is the same as that derived by 
iGriffiths and Tavara 1)19981 ) For selected a lleles, the frequency spectrum cannot be 
obtained by the method of lGriffitha l|2003l) except in the special case in which the 
selection intensity is inversely proportional to the population size at all times in 
the past. 

The freguencv spectr um of alleles closely linked to selected loci is also of in- 
terest. iBraverman et al. ( 1995 ), Fav et al. ( 2002 ), Kim and Stcphan (2002), and 
others have simulated the effects of selected sites on the frequency spectrum of 
closely linked neutral sites, with the goal of finding evidence of background se- 
lection against deleterious mutations and genetic hitchhiking caused by positive 
se lection of advantageous m ut at ions . 

IWilliamson et alJ pOOR^ recently considered the combined effects of population 



growth and selection on the frequency spectrum. Their model was of a popula- 
tion that was of a constant size until r generations in the past, at which time it 
grew instantaneously by a factor v and remained at the new size until the present. 
IWilliamson et al. (2005) developed a likelihood method for estimating both t and 
u from the spectrum of sites assumed to be neutral and for estimating r, i' and 
a scaled selection intensity 7 for non-neutral sites. Their method was based on 
numerical solutions for the frequency spectrum using both Kimura ( 1955) series 
solution for neutral alleles and numerical solutions to the backwards equation for 



NON-EQUILIBRIUM FREQUENCY SPECTRUM 



selected alleles. IWilliamson et alJ 1 20051) applied their method to a previously pub- 
lished data set of 301 genes in the human genome and found evidence both of 
population gro wth and of purify ing selection at non-synonymous sites. In a related 
study, Bustamante et alJ l|2005|) found evidence of differences in selection intensity 
among different classes of genes in a data set for more than 6000 loci in humans. 

In this paper, we will explore in more detail the allele frequ ency spectru i n in a 
population of variable size. Our goal is similar to that of Williamson et al.l 1(20051) 
in modeling the combined effects of selection and population growth. We derive the 
frequency spectrum from the forward equation, first for a Markov chain and then for 
a diffusion approximation to that Markov chain. The forward equation provides a 
natural way to compute the spectrum as it approaches an equilibrium from an arbi- 
trary initial condition and to model the time-dependence of the spectrum resulting 
from the time-dependence of population size and selection intensity. Furthermore, 
the forward equation provides a way to incorporate the effects of immigration. We 
show that our approach recovers known results for an equilibrium population and 
present some numerical results for an idealized model of recent human populations. 

2. Markov chain 

The model is of a monoecious randomly-mating diploid population containing 
N(t) individuals at time t, which in this section takes integer values representing 
discrete non-overlapping generations. We assume a large number of identical and 
independent loci. At each locus there are only two alleles A, the derived allele, and 
a, the ancestral allele. In generation t, the set of loci is described by the row vector 
with j*^ element fj (t) that is the expected numbers of loci at which A is found on 
j chromosomes, 1 < j < 2N{t). Thus f2N{t){t) is the expected number of loci fixed 
for A in generation t. The model assumes that the pool of loci fixed for a is so 
large that it can be assumed to be not reduced b y the creation o f polymorphic loci 
by mutation - the infinitely many sites model of lKimural i|l969l) . 

The change in fj {t) because of genetic drift and mutation is described by the set 
of difference equations 

2N{t) 

(1) f,{t + 1) - ^ Mt)P^At) + Mj{t), 1 < J < 2N{t + 1). 

1=1 

The first term on the right hand side represents the combined effect of genet ic drift 
and na tural selection on loci that are already polymorphic: in the notation of lEwenal 
l|2004l) . Pij{t) is the probability that a locus with i copies of A in generation t will 
have j copies in generation i + 1. The pij (t) are easily derived for the Wright-Fisher 
and other models l|Ewensl IJOO^ . 

The second term on the right hand side represents the creation of new poly- 
morphic loci by mutation and immigration. The influx of mutations is modeled by 
assuming that each of the 2N{t) copies of a at a monomorphic locus mutates to an 
A with probability /x per generation. Therefore 

(2) Mj{t)^2N{t)^i5i^j 

for mutation alone, where (5i.j = 1 if j = 1 and otherwise. Under the infinitely 
many sites model, the mutation rate is assumed to be so low that the effect of 
mutation on loci already polymorphic can be ignored. However, mutation can be 
incorporated into pij (t) if necessary. 
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Immigration from another population can be accounted for both by modifying 
Pij to aUow for the effect of immigration on loci that are already polymorphic 
and Mj{t) to allow for the creation of new polymorphic loci. Immigration, unlike 
mutation, can create polymorphic loci for which j > 1 immediately. In this paper, 
we will restrict our analysis to the case with mutation only. 

Given an initial frequency spectrum, /j(0). Equation Q can be iterated to obtain 
the frequency spectrum at any time in the future. Even if N is constant, there is 
no equilibrium solution to Equation (Q because the number of loci fixed for A will 
increase each generation. However, an equilibrium solution for /j, I < i < 2N — 1, 
the spectrum for polymorphic loci, does exist and is found by solving the matrix 
equation 

(3) f = fP+{^ 

where P is the {2N — 1) x {2N — 1) matrix with elements pij for 1 < i, j < 2N, e is 
a row vector with the first element 1 and the remaining elements 0, and 9 ~ AN/i. 
The solution is 

(4) f=leiI-P)~\ 

where / is the (27V — 1) x (2iV — 1) identity matrix. 

This result is equivalent to that obtained using the backwards equation for the 
Markov chain. The frequency spectrum of polymorphic loci is proportional to the 
sojourn times (iij in the notation o f Ie wend 1120 04 ll. Equation Q) is obtained from 
the solution to Equation (2.143) of Ie wens! 1)20041 by multiplying by 9/2, which is 
the rate of influx of mutations each generation. 

The advantage of the formulation presented here is that it also describes the 
approach to the equilibrium. The rate of approach depen ds on the s econd largest 
eigenvalue of P, which for a neutral allele is 1 ~ 2^ (,Ewensl l2004f) . The mean 
number of loci fixed for A increases by 

2N-1 

(5) E /^(^>^.2jv 

per generation, which reduces to 2Nfj,Pi at equilibrium, where Pi is the probability 
of fixation of each mutant. 

3. Diffusion approximation and new boundary condition 

Let us start by considering the time-homogeneous case with no mutation from 
the ancestral type, but where we can start at time with some derived alleles 
already present. Because we want to eventually allow varying population sizes, 
assume that that the population is described by a time-homogeneous Markov chain 
with state-space {0, 1, ... , 2Np}. Suppose that if we shrink space by a factor of 
2Np and speed time up by a factor of 2N, then this chain converges to a diffusion 
process on [0, 1] with generator Q = a{x)-j^ + ^b{x)-^ for appropriate coefficients 
(this is the scaling regime that is appropriate for models such as the Wright-Fisher 
chain with or without selection). 

Suppose at time that there are countably many loci at which derived alleles are 
present, with respective (non-random) frequencies xi, a;2, . . .. Once we have passed 
to the diffusion limit, the frequency spectrum at time t is just the intensity measure 
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(that is, the expectation measure) of the point process that comes from starting 
independent copies of the diffusion process at each of the Xi and letting them run 
to time t. In other words, the intensity measure is obtained by taking the sum 
of point masses and moving it forwards an amount of time t using the semigroup 
associated with the generator Q. 

More generaUy, if the initial configuration of frequencies is random (so that it 
can be thought of as a point process on (0, 1)), then the frequency spectrum at time 
t is obtained by taking the measure that is the intensity of that point process and 
again moving it forwards an amount of time t using this semigroup. 

For i > the resulting measure will be absolutely continuous with respect to 
Lebesgue measure and have a density f°{y,t) at frequency y e (0, 1). We will also 
refer to this density as the frequency spectrum. It is immediate that /° satisfies the 
Kolmogorov forward equation equations that go with the generator Q (with initial 
conditions corresponding to the intensity measure of the point process of initial 
frequencies). That is, 

(6) |r(2/,i) = -^[a{y)ny,t)] + l£^[b{y)riy,t)], 

with limyj^o /"{y, t) and limj,|i f°{y, t) both finite and appropriate boundary condi- 
tions at i = (in particular, if the point process of initial frequencies has intensity 
h{y) dy, then limtxo /°(2/, t) = h{y)). 

Now we want to introduce mutation from the ancestral type as time progresses. 
In the Markov chain model, this corresponds to new mutants arising in the popu- 
lation at rate ^P P^r unit of Markov chain time, where 9 is independent of N. The 
initial number of mutants at a locus is 1. This is equivalent to mutants appearing 
at rate 27V|jO per unit of rescaled time, with the initial proportion of mutants at a 
locus being ^. 

Imagine now that we pass to the diffusion limit for the allele frequencies, but 
for the moment still work with a finite N for the description of the appearance of 
new mutants. That is, we think of our evolving point process as having new points 
added at location 2^ at rate ^2Np, and that the locations of these points then 
evolve as independent copies of the diffusion with generator Q. 

We will make substantial use of t he theory of entrance la ws for one-dimensional 
diffusions laid out in Section 3 of IPitman and Yoil l)l982J) . Write Pt{x,dy) for 
the semigroup associated with^. Th is is the semigroup of the 0-diffusion in the 
terminology of lPitman and Yon 1 1982|) . The contribution to the frequency spectrum 
from mutations that appear after time is 

(7) '4'^ Ij^- {2k' '')''■ 

Choose a scale function s for the 0- diffusion such that s(0) = (so that s is then 
unique up to a positive multiple). As lPitman and Yon l)l982r) observe, 

(8) P^{x,dy):^-^Pu{x,dy)s{y), 0<x,y<l, 

is the semigroup of a diffusion that never hits (this 1 -diffusion is the Doob h- 
transform that corresponds to the naive idea of conditioning the 0-diffusion never 
to hit 0). Moreover, this semigroup can be extended to allow starting at by setting 
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(9) P^{0,dy)=lhnPlix,dy). 



The resulting extended process can start at but it will never return to 0. 

Assume now that s'(0) > 0, which will be the case in the diffusions that are of 
interest to us. We can choose the free multiplicative constant in the definition of 
the scale function s so that lim^io s{y)/y = s'(0) = 1. Then 

(10) ^in.JN,P. (^,.,) . 5liM^ =., A„(,,) 



in the notation of Pitman and Yor. Thus 

(11) ;i^JN{p[Pt-s [^^.dy) ds ^'-j\,^.{dy)ds .: a>.(d,), 

say. 

As observed in lPitman and Yon l|1982|) . the family (A„)„>o is an entrance law for 
the semigroup of the 0-diffusion, and so it has densities that satisfy the Kolmogorov 
forward equation associated with the generator Q — intuitively, (A„)„>o describes 
that injection of an infinite amount of mass at location at time 0, with this 
mass subsequently evolving in (0, 1) according to the dynamics of the 0-diffusion. 
Consequently, the family ($t)t>o also satisfies the Kolmogorov forward equation 
associated with the generator Q — again intuitively, (<l>t)t>o describes a continuous- 
in-time injection of mass at location 0, with this mass again subsequently evolving 
in (0, 1) according to the dynamics of the 0-diffusion. That is, if we write (pt for the 
density of $t, we have that 

(12) |<^*(2/) = -^Hy)<t>t{y)] + \^Hy)Uy)]- 

It remains to work out what the boundary conditions for (pt are. Following 
IPitman and Yon l|1982r) . introduce the [-diffusion^ which is the 0-diffusion condi- 
tioned to hit before 1. The J,-diffusion has the Doob /i-transform semigroup 

(13) p/(^,dy)=fi_4^') ' P,{:,,dy)ll-'^y^ 



s{l)J - ' -V s{l) 

From IWilliamj l)1974f) . the t-diffusion started at and killed at the last time it 
visits 2/ > is the time-reversal of the |-diffusion started at y and killed when it 
first hits 0. Write {Qi)t>o for the semigroup of this killed process. Since we have 
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normalized the scale function s so that s{y) « y for y close to 0, 
\im y(t)t{y) = lim s{y)(j)t{y) 

yiO yiO 



yio ^^'Jn 2 dy 



6 f'PlMdy), 

= — iini / as 

2 vio Jo dy 

2 vlo Jo dy 

2 !/io Jq dy 

2 vlo Jo dy 

Note that if {Bt)t>o is a standard Brownia n motion and T :— inf{i > : i?* = 0}, 
then, by Equation (3.2.1) and Section 3.1 of|Knightl fl981), and 

\-KS 

2A ^2a'''""-^'"' 



dy 
(15) 




= 2y. 
Observe also that a scale function for the J,-diffusion is 

,^„s / N •s(x)s(l) 

(^^) ^(^) ^- ^(TFTR- 

By standard one-dimensional diffusion theory, if we compose the killed J,-diffusion 
with (T, then the resulting process is a time-change of standard Brownian motion 
killed when it first hits 0, with the time-change given by the c orresponding speed 
measure (see, for example, V.7 of iR.ogers and Williamsl lj200(T 'l. Moreover, since 
(t(x) ~ a; for x close to 0, the speed measure to for the |-diffusion satisfies m{dx) ^ 
dx/b{x) for X close to (beware that the definitions of the speed measure can vary 
from author to author by m ultiplicative constants, we are using the definition of 
iRogers and WilUameJ l|2nn(t ). Therefore 

(17) ^imr%M,,4n^^.,H^_^. 
2 aio Jo dy 2 yio b{y) yio h{y) 

Write f{x,t) for the frequency spectrum of the model with mutation from ances- 
tral type. We have f{x, t) — f°{x, t) + (l)t{x), where /° is defined for the appropriate 
initial conditions at t = 0. If we want to start with all alleles ancestral type, then 
the initial conditions at t = are null and /° = 0. Combining what we have 
obtained above, we find that 

(18) p{x,t) = ~-^[a{x)f{x,t)] + ^^[b{x)f{x,t)]. 

with appropriate boundary conditions at t = (in particular, \imtio f{x,t) = 
if we start will all alleles being ancestral), and further boundary conditions 
lim^ixo xf{x, t) = 6'lim2;xo -^ and Yivcixn f{x, t) finite. 
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Now consider a timc-inhomogeneous diffusion with generator a{x, t)-j^ + '^b{x, t) -j^ 
and suppose also that p is now a function p{t) of time. By first considering the case 
where a, b and p are piecewise constant, using the above analysis, and then taking 
limits, we get that the frequency spectrum solves 

(19) |/(y,0 = -^Ky,i)/(y,t)] + \^^[b{y,t)f{y,t)] 

with appropriate boundary conditions at i = and further boundary conditions 
limaio2//(2/,i) = OVnn^^^a ^^ and lim^^ii f{y,t) finite. 

For the purposes of a numerical solution, it is more convenient to consider the 
function g{x,t) :— x{l — x)f{x,t) which satisfies 



(20) lg{:c,t)^-x{l-x)^ 



a{x,t) 



x{l — x) 



g{x,t) 



2 



x{l — x) d 
2 fe2 



b{x,t) 
x{l — X 



g{x, t) 



with appropriate boundary conditions at t = and further boundary conditions 
lim^XO 9{x, t)^e lim^io -g^ and lim^^ i g{x, t) == 0. 

As an example, consider the case where a{x) = Sx{l — x), and b{x) — x{l — 
x)/p{t). This is a Wright-Fisher diffusion with selection and varying population 
size. The corresponding forward equation is 

(21) |.9(., t) = -Sxil - x)l [,(., .)] + :^^ [,(., t)] 
with boundary conditions 

(22) lim5(a;,t) == dlim^^P^ = 0p(t). 
^ ' xiO^^ ' ^ xio x{l -x) ^^ ' 

4. Equilibrium solution 

When p(t) = 1 is a constant and the coefficients a and b do not depend on time, 
then we expect f{-,t) to converge to an equilibrium / as i — > oo. 
From Equation (|19|l . the function / should satisfy 

(23) = -^[a{x)fix)] + l^[bix)fix)] 

with boundary conditions lim2:xo xf{x) — lima;|o j^ and lima;|i f{x) finite. 

In order to solve this equation, suppose first that the diffusion process is on 
natural scale in [0, 1], so that a = 0. In that case. Equation (|23|l becomes 

(24) 0=i^[6(x)/»], 

so that b{x)f{x) — cq + cix for some constants cq and ci. Assume that < 
lim:r|i(l — x)/b{x) < oo. Satisfying the boundary conditions requires that f{x) = 
e{l-x)lb{x). 

Suppose now that the diffusion process is not on natural scale and that s is a 
scale function with s(0) = and s(l) ~ 1, so that the image of the diffusion under 
s is a diffusion on [0, 1] in natural scale. The generator of the image diffusion is 

(25) \[{s'os-^){x)f{bos-^){x)^, 
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and hence the image diffusion has the associated frequency spectrum 

^^^^ ^^^^ ■" [(s'os-i)(a;)]2(&os-i)(a;) 

from what we have just observed. It follows from the usual change of variable for- 
mula for densities that the original diffusion has the associated frequency spectrum 

(27) fix) = ihos)ix)s'{x) = ^-^^^^^- 
In terms of coefficients, 

(28) 5(^) = r exp (-2 J^ a{z)/b{z) dz) dy 

Jo exp (-2 J^ a{z)/b{z) dz) dy. 
Also, 

1 dm, , 

s'(xjo(x) rfa; 

where m is the speed measure corresponding to the scale measure s (recall that w e 
are using the normalization of the speed measure in lRosrers and Williama Ij200(l )). 
Thus 

(30) f(:c)=ePo{x) — {x), 

where Po{x) is the probability the diffusion will be absorbed at given that it 
starts at x. This well-known equ ation can be es tablished directly via an ergodic 
argument using time-reversal - see lGrifhthsl l|2003|) for a discussion of the history of 
this technique. Note that our derivation of Equation 119() also used time-reversal. 
For example, when a = and b{x) = x{l — x) (the neutral Wright-Fisher diffu- 
sion), /(x) — 9/x, agreeing with Equation (9.18) of Ie wend (|200J). Similarly, when 
a{x) = Sx{l — x) and b{x) = x{l — x) (the Wright-Fisher diffusion with selection 
and no dominance). 



e 



2S h _ „-2S{l-x) 



(1 



^^'^ ^^")-" (e2^-l)x(l-x) ' 

agreeing with Equation (9.23) of lEwensI 1 20041) . 

5. A SYSTEM OF ODES FOR THE MOMENTS IN A WRIGHT-FISHER DIFFUSION 
WITH VARYING POPULATION SIZE 

Suppose in this section that a{x) — Sx{l — x) and b{x) = x{l ~ x)/ p{t). This 
is a Wright-Fisher diffusion model with time-varying population size and constant 
selection and without dominance. Equation H21|l applies. 

Put /i„(i) := /q x"g{x, t) dx for n = 0, 1, 2, . . .. Integrating by parts, we get 

"1 d 

x"a;(l — x)-—g{x, t) dx 
c/a; 

(32) =[(x"+i-."+2).9(.,0]; 

- / {{n+l)x'' -{n + 2)x''+^)g{x,t)dx 
Jo 
= {n + 1)^„ - {n + 2)^„+i(t) 
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Similarly, 

x"'x{l — x)—-^g{x,t) dx 

(x"+i-x"+2)|-5(x,t)l 
dx Jo 

^ d 

((n + l)x" - (n + 2)a;"+i)— .9(2;, i) dx 
ox 

((n + l)x" - (n + 2)a;"+i)— .g(a;, t) dx 



(33) 



9a:; 

= -[((n + l)a;"-(n + 2K+i)5(x,t)]; 

+ / ((n+l)nx""H{n7^0}- (n + 2)(n+l)a;").g(a;,t)dx 
Jo 
= [l{n = 0}^p(t)] + [(n + l)n^lr^-lit)l{n ^ 0} - {n + 2){n + l)fin{t)] 

We thus get the coupled system of ODEs 
(34) Mo (^) = ^ - ^Mo (t) + 5 (mo (i) - 2mi (t)) 

and 

mUO = ^^ [(^ + l)n^in-l{t) - (n + 2)(n + l)M„(t)] 

+ S{{n+ l)M„(t) - (n + 2)M„+i(t)) , n > 1. 

When 5' = (that is, the neutral case) this system is lower triangular and we 
can be solved explicitly. The ODE for /io has solution 



(36) Mo(t) = Mo(0) exp ( - / — ds ^ ' " ^' "'' ^^° " '"^ '^"^ '^^ 



/o pW ; 2 exp(/„*^d. 
Given Mn-i, the ODE for fin has solution 

M„(i)=Mn(0)cxp(^-(^" + ^)^ -i^rf. 



^'') /n + 1\ /o ^A^"-i(-) -P ((-r) /o^ ^ du) 



ds 



-p((T)/o;;fe'^^) 



We can draw some conclusions from these equations about the effect of p on the 
asymptotic behavior of /in- For example, recall that the measure f{x,t) dx is the 
intensity of the point process on (0, 1) that records the frequencies of derived alleles 
at all the loci at which derived alleles are present with non-zero frequency at time 
t. Hence 2/io(i) = 2/q x{l — x)f{x,t)dx is the expected value of the total of the 
heterozygosities summed over all loci. For any initial conditions and any p such 
that /„ —7-p:d,t < 00, the expected total heterozygosity 2/io(i) is asymptotically 
equivalent to 9t as t ^> 00. 
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6. Explicit recurrences for the moments in a Wright-Fisher 
diffusion with exponentially increasing population size 

Suppose in this section that p{t) = e^* with R > 0, a{x) — 0, and b{x) — 
x{l—x)/e^^. This is a neutral Wright-Fisher diffusion model with constant selection 
and exponentially increasing population size. We will obtain an explicit recursive 
recipe for the moments iJ,n(t). For simplicity, suppose that f{x,0) = 0, so that 
/i„ = for all n. A similar development holds for other initial conditions. 

Recall that the exponential integral function Ei is given by Ei (x) = — J_ i~^e~* di, 
where the principal value is taken if x > (although we are only interested in the 
case X < 0). For a; < 0, Ei (x) = — r(0; —x), where F is the usual upper incomplete 
gamma function. For x > 0, 

r e"* - 1 
Ei {-x) = 7 + log(2;) + / dt 

Jo i 

^^^^ ^ {-xV 

= 7 + log(x)+}_^^-^, 

where 7 is Euler's constant l|Gradshtevn and RvzhikL |200(]|) . 
For n = 0, 

B e-«* / / 1 
(39) ^o(i) = TTTje^^ Ei -- ) - Ei 



2R \ \ R) \ R 

For n = 1, 2, . . . define a linear operator $„ by 

so that fin = ^nfJ.n-1 = ■ ■■ = <l>„$„_i • • • $1^0- Set 



hk{t) = exp 



R 



(41) r . n ..... . .,_ , .s _-H, 



X 



Ei(-('^+2UUEir-r+'^^" 



2 J RJ \ \ 2 J R 

so that /io = 2i7^o- It follows from a straightforward integration that 

(42) ^nhkit) = -^^^^ll^[h,it) - h,M- 

Hence 

n 

(43) Unjt) := y^^Cn^khkjt), 

k=0 

where co,o — td ^^'^ ^^^ other Cn,k are given recursively by 

(44) Cn,k = 7;r+2y37fc+2yCn-i,fe, l<i<n-l, 
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and 

(45) 

For example, 
(46) 
and 

(47) 



n-l 

' 2^ 'fn+2 
k=0 
n-l 



CV) 



Z^ ('n+2\ _ (k+2\ 

k=o V 2 ; V 2 ; 

n-l 



Cn-l.fe 



k=0 

Clfi = 
Cl ,1 = 



1 1 

3-12^ 
-1 

4R 



1 



1 



C2,0 



C2.1 



6- 14i? 

3 (-1) 
6-3 4i? 



3 
20R 



1 

Ir 






1 

Ir 



lOi? 



Set ^(x) := E°li Tn- Then 



/i/j (f ) = exp 



(48) 



Ei 



fc + 2 
2 
k 



exp 



fc + 2 
2 



R 
2 ji? 



Ei 



fc + 2' 



-flt 



i? 



X Rt + tp 

It follows from this that 
(49) 



R 

fc + 2 



A*o(t) 



V' 



fc + 2\ e 
2 



-_Rt 



i? 



0t 



when R is large (recall that 2y^o(i) is the expected total heterozygosity summed over 
all loci). For n > 1, the two observations that exp ( — C^j^) ^-^ 1 will be very close 
to 1 for even moderate values of R and that ^j, Cn,k — show that the contribution 
to fin (t) from the Rt term will almost cancel out and the primary contribution will 
be from the ip terms. 

7. Frequency spectrum in a finite sample 

The function f{x,t) approximates the frequency spectrum in a very large pop- 
ulation. In a sample of n chromosomes, we can observe only the finite spectrum, 
fi{t), which is the distribution of the number of chromosomes at which there are i 
derived alleles {0 < i < n). In this context, fi{t) is similar to fi{t) defined for the 
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Markov chain formulation, but here t is continuous. The finite spectrum is obtained 
from f{x,t) by assuming samphng with replacement at each locus independently; 

(50) m) = (A ^ x\i- xr-'f{x, t) dx 

iGriffithal pOO,*^ . We can express this equation in terms of the moments of g{x, t) = 
a;(l — x)f{x, t) about x = 0: 

(51) /,(t)= (^';)"g\-l)^(^"~^^"^)^,+,_i(t) 

This expression involves an alternating sum and so, from a numerical point of 
view, it might be preferable to have a system of ODEs for the fi themselves rather 
than for the Hi. Unfortunately, a system of ODEs for the fi appears to be rather 
complicated: it is doubly indexed by i and the suppressed index n and, moreover, 
the alternation present in Equation (|51l) is just "transferred" to the coefficients of 
the ODEs. 

For use later in Section|^we make the following small observation. Suppose that 
x I— » f{x,t) is decreasing. Observe for 1 < i < n — 1 that, by an integration by 
parts, 

(52) 

/,;+l(t)-/,(i) 



1 

i+l 



(t + iy.in-iy.j. 



x'+^{n - i){l - x^-'-^ - (i + l)x'(l - xY-'\ f{x, t) dx 



— [x'+\l^xr~']f{x,t)dx 



{i + l)\{n-i)\ Jq dx 

' x'+\\-xY~'—f{x,t)dx 



{i + l)\{n-i)\ J^ ' ' dx 

<0 

(there are no boundary terms because of the boundary conditions on f{-,t)). Thus 
fi{t) is decreasing in i when x i— > f(x,t) is decreasing. 

8. Numerical analysis 

If the population size or selection intensity vary only somewhat with time, nu- 
merical solutions to Equation 1)21(1 can be obtained with standard methods. The 
function NDSolve in Mathematica (version 5) and the function pdepe of MATLAB 
(version 7) both provide solutions using default settings for those programs. With 
extreme population growth, as in the model of human history we consider in Sec- 
tion O neither of these programs provides accurate solutions, so it was necessary 
to write a program tailored to the problem. 

Large gradients at the boundary a; = make it necessary to introduce a non- 
uniform grid with A'^ internal points. With the view toward integrating the numeri- 
cal solution in order to compute moments and the finite spectrum we select the grid 
with the smallest spatial increment Aq = {xi — 0) of the order 10~^ and the largest 
Aat = (1 — xn) ^ 10^"^. The spacing is kept constant for a few nodes Xi, ...^Xi+k 
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and then doubled so that 

Xi+k+i — Xi+k = A,;+fe = 2Ai+fc_i = 2{xi+k — Xi+k-i). 

This process is repeated until the spacing reaches the maximum size An ^ 10"^. 
The numerical domain is thus separated into sub-domains within each of which 
the spacing is uniform. This guarantees that centered difference schemes used in 
the sub-domains give a second order accurate approximations to the corresponding 
differential operators. For the right end-point of the uniform sub-domains, x^+fc, 
(referred to as the edge-node) the values at Xi+k~2 and Xi^^+i can be used, since 

Xi+k — XiJ^k-2 = Xi+k+1 — Xi+k- 

li t — Atfc and x = Xi we denote the numerical approximation to g{x,t) by G^. 
A standard centered second order differencing is used to approximate the diffusive 
term and one-sided first order difference with the direction depending on the sign 
of S (so-called upwinding) is used for the advective term. For the time-stepping 
an implicit Backward-Euler scheme with a fixed step-size At — An is used. The 
following system of algebraic equations is obtained: 

/ a;j(l - a 



(53) G^+^ = G'l -f At 



\2p{t){A,] 



Gk-\-\ C\^ik+1 I y-^fc + l 

i-Edge ~ '^^i + ^i+l 



_oA I Xi(l Xj) v„k+UwR _ ^k+UwLl 

where Edge — 2 if Xi is and edge-node and Edge = 1 otherwise; UwR — 0,UwL — 
— 1 if 5* > and UwR — 1, UwL = if 5 < 0. The boundary conditions are 
Gg = giAt * k) and G%^, = 0. 

Overall the truncation error is of the order Ajv- In each time step, in order to 
solve the system (|53l ') it is necessary to invert a sparse matrix with a large condition 
number {cN « 10*), hence a gain in accuracy that would come from decreasing the 
time step is offset by the loss of precision in the inversion of such a matrix. 

9. Model of recent human population growth 

We considered a highly simplifie d model of the histo r y of p opul ation sizes of mod- 
ern hu mans similar to that used bv lReich and Landed l|200l|) and lWilliamson et al.l 
l|2005|) but not requiring the assumption of an instantaneous change in popula- 
tion size. We assumed a stable population containing A^o = 10, 000 individuals 
until 150, 000 years ago, which we took as time i = 0. We assumed a genera- 
tion time of 25 years. We measured time in units of 2A*'o, so the present is at 
t = 6000/20000 = 0.3. At i = the population began to increase in size exponen- 
tially at a scaled rate R = 2Nor, where r is the exponential rate per generation. 
We assumed additive selection with heterozygous fitness 1 + s relative to aa ho- 
mozygotes, with the selection coefhcient is also scaled hy Nq: S = 2Nqs. The value 
R — 40 corresponds to a current size of 1.63 x 10^. Reich and Lander assumed a 
current size of 6 x 10^, but that did not take account of the fact that the effective 
size of human populations is roughly 1/3 of the census size lHil]L ll972j) . 

We assume that the spectrum at t = is the equilibrium spectrum: 

(54) g{x,0) = 9il~x) 
for 5* = and 

(55) g{x, 0) = 9 ^-5^- ^ 
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otherwise. The numerical solutions for f{x, t) at times t — Q and t — 0.3 are plotted 
in Figures 1-3 for the respective choices 0, +2, —2 of the selection parameter S (the 
mutation parameter is taken to be = 1 - a different choice of 9 merely rescales 
the spectrum). 

FIGURE 1 HERE 

FIGURE 2 HERE 

FIGURE 3 HERE 
For neutral alleles, the equation for the 0th moment of g{x, t), which is half the 
expected total heterozygozity, can be solved exactly. It is of interest to separate 
the solution into two parts, one representing alleles present at t = (old alleles) 
and the other representing alleles that arose by mutation after t = (new alleles). 
For old alleles, 

(56) — -— = -e ^ fio.o 

at 

with initial condition ^0,0 (0) = | and for new alleles 

(57) -^-2-' ^°- 

with initial condition /io,n(0) = 0. Equations (|5(j|l and (|57|l can be solved in 
terms of exponential integrals. With R — 40, we find /io,o(0.3) = 0.496' and 
Mo,n(0.3) = 0.156*, which implies that under this model, roughly 30% of the ex- 
pected heterozygosity at neutral sites is attributable to mutations that arose in the 
past 150,000 years. 

For selected alleles (5' 7^ 0), the system for the moments is not closed. An ap- 
proximation to the moments can be obtained by truncating the system and setting 
the first neglected term to it's initial value. Alternatively, the moments can be 
computed by first solving for g{x,t) and numerically integrating. Both approaches 
present certain numerical difficulties. The plots in Figure 4-6 were obtained for a 
sample of size n = 20 by numerically solving the truncated system with 160 equa- 
tions using MATLAffs ode45 (for S* = 0) and odelBs (for S = ±2) routines. For 
the ne utral case 5* = there is a simulation algorithm due to lGriffiths and Tavara 
lll998|) for approximating the finite spectrum, and we compare that approximation 
with the results from the numerical solution of the the system of ODEs in Figure 
7 for a sample of size n = 40. The mutation parameter for Figures 4-7 is taken 
to be 6* = 1. Again, a different choice of 9 merely rescales the finite spectrum. It 
appears from Figure 1 that the frequency spectrum at time i = 0.3 is a decreasing 
function and hence, from the remark at the end of Section [7| the corresponding 
finite spectrum plotted in Figure 7 should also be decreasing. The undulations in 
the plot are therefore numerical artifacts. 

FIGURE 4 HERE 

FIGURE 5 HERE 

FIGURE 6 HERE 

FIGURE 7 HERE 
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Figure 1. Frequency spectrum f{x,t) ~ g{x,t)/{x{l — x)) at 
times t — and t = 0.3 with parameter values i? = 40 and S = 0. 
Obtained by numerically integrating the PDF. The values of / are 
restricted to the interval [0, 100]. 
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Figure 2. Frequency spectrum f{x,t) = g{x,t)/{x{l — x)) at 
times t — and t = 0.3 with parameter values R — 40 and S — +2. 
Obtained by numerically integrating the PDF. The values of / are 
restricted to the interval [0, 100]. 
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Figure 3. Frequency f{x,t) — g{x,t)/(x[l — x)) at times t — Q 
and t = 0.3 with parameter values i? = 40 and S = —2; Obtained 
by numerically integrating the PDF. The values of / are restricted 
to the interval [0,100]. 
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Figure 4. Finite spectrum for a sample size n = 20 at times t — 
and t — 0.3 with parameter values i? = 40 and 5* = 0. Obtained 
by numerically integrating the system of ODEs for the moments 
oig{;t). 
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Figure 5 . Finite spectrum for a sample size n = 20 at times t — 
and t — 0.3 with parameter values R — AO and S = +2. Obtained 
by numerically integrating the system of ODEs for the moments 
oig{;t). 



24 STEVEN N. EVANS \ YELENA SHVETS ^, AND MONTGOMERY SLATKIN ■ 




f.{n=20,t=0) 
* f.{n=20,t=0.3) 



0.3 0.4 0.5 0.6 0.7 

i/n = number of derived alleles/ sample size 



0.8 



0.9 



Figure 6. Finite spectrum for a sample size n = 20 at times t — 
and t = 0.3 with parameter values i? = 40 and S = —2. Obtained 
by numerically integrating the system of ODEs for the moments 
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f(n=40, t=0.3) from simulation aigorithm 
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Figure 7. Finite spectrum for a sample size n = 40 at times t — 
and t = 0.3 with parameter values i? = 40 and S ^ Q. Obtained 
by numerically integrating the system of ODEs for the moments of 
g{-,t). Compared at time t = 0.3 with results from the simulation 
algorithm of Griffiths and Tavare. 



